%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O-Ring Production Networks
%   Author: Demir, Fieler, Xu, Yang
%   Last updated: January 2021
%
% Purpose: Endogenous Targeting (step 3)
%   - Select two models to report model fit
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; diary off; clc; clear; 

    % load estimates (exogenous targeting)
    load('../../output/estimates/fullmodel_estimates.mat', 'est');            
 
    % load endogenous targeting calibration summary
    load('../../output/temp/endogtarget_initeqm_grid.mat')     
   


%% Select two models
for nc = [1 17]
tic    
    
    %% Pick best fit calibration
    
        % best fit nuv index
        fval_matrix_adj = fval_matrix + 100*(converg.iter_outer==100*ones(Nc,Nv)) + 100*(converg.iter_inner==500*ones(Nc,Nv));
        [~, nv] = min(fval_matrix_adj(nc,:));
                
        % nu_c and nu_v
        nu_c = nuc_matrix(nc, nv);
        nu_v = nuv_matrix(nc, nv);
        fprintf('\nCounterfactual(%.0f,%.0f): nu_c = %.4f, nu_v = %.4f\n', nc, nv, nu_c, nu_v)
        fprintf('Saved calibration: fval = %.4f, err_outer = %e, err_inner = %e, iter_outer = %.0f, iter_inner = %.0f\n', ...
           fval_matrix(nc,nv), converg.err_outer(nc,nv), converg.err_inner(nc,nv), converg.iter_outer(nc,nv), converg.iter_inner(nc,nv));    
        
        
    %% Initial equilibrium
    fprintf('\nInitial Equilibrium\n');
    
        % parameters
        define_param
        define_param_est_endogtarget     
        param.nu_c = nu_c; %calibrate nu_c
        param.nu_v = nu_v; %calibrate nu_v

        % data moment
        define_data_moment

        % constants
        define_omega_dist 
        define_const_endogtarget 

        % solve initial equilibrium
        eqm = solve_eqm_outer_endogtarget(param, const);    
        define_wage_schedule
        const.phi_v = eqm.phi_v; %add phi_v to const

        % calculate model moments and objfcn
        model = cal_model_moment(param, const, eqm); 
        bartik = cal_chg_bartik(param, const, eqm);     
        fval = cal_objfcn(data, model, bartik);

    
    %% Store results
    
        % display
        fprintf('\nparam.nu_c = %.2f, param.nu_v = %.4f, fval = %.4f\nInitial equilibrium: n(qmax) = %e, err_outer = %e, err_inner = %e, iter_outer = %.0f, iter_inner = %.0f\n', ...
            param.nu_c, param.nu_v, fval, eqm.nq(param.Q), eqm.err_outer, eqm.err_inner, eqm.iter_outer, eqm.iter_inner);     
    
        % store
        result.param = param;
        result.const = const; 
        result.eqm = eqm;
        result.model = model;
        result.bartik = bartik;
        
        % save
        if nc == 1
            result1 = result;
        elseif nc == 17
            result2 = result;
        end

toc
end



%% Output tables
diary off

    % create folder
    mkdir('../../output/tables')    


    % endogenous targeting model fit
    save('../../output/temp/endogtarget_select_modelfit.mat', 'result1', 'result2');

    % exogenous targeting model fit
    run('../1 estimation/main_fullmodel_moments.m')
    clearvars -except data param model bartik

    % 
    load('../../output/temp/endogtarget_select_modelfit.mat', 'result1', 'result2');
    delete('../../output/temp/endogtarget_select_modelfit.mat', 'result1', 'result2');


    % output
    out_endogtarget_tables




 
    